library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
## 
##     expand
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loaded glmnet 2.0-18

Note on how data are structured in each csv: First 30 are non drug Next 30 are drug

Alternate between left and right turning per row

Please be sure to put the data files under the /data folder

### Read in all the datafiles
r1<-(read.csv('data/r1.csv', header = FALSE))
r2<-(read.csv('data/r2.csv', header = FALSE))
r3<-(read.csv('data/r3.csv', header = FALSE))
r4<-(read.csv('data/r4.csv', header = FALSE))
r5<-(read.csv('data/r5.csv', header = FALSE))
r6<-(read.csv('data/r6.csv', header = FALSE))
r7<-(read.csv('data/r7.csv', header = FALSE))
r8<-(read.csv('data/r8.csv', header = FALSE))

## Create two vectors for outcome variables based on how we know the data are structured
drug_cond<-c(rep(0,30),rep(1,30))
dir<- c(rep(c(0,1),30))
y_colnames = c('drug_cond','dir')
### Function returns mean cross validated accuracy for glmnet models for each rat
### So 8 rats * 6 cv splits* 2 accuracy scores (drug condition/direction)

manual_cv<- function(df){
running_acc_dg = {}
running_acc_dir = {}
for (i in 0:5) {
  df=cbind(df,drug_cond,dir)
  s=(10*i)+1
  e=10*(i+1)
  test_set = c(s:e)
  fit_df = df[!(rownames(df)%in%test_set),]
  test_df = df[test_set,]
  fit_dg = fit_df[,'drug_cond']
  fit_dir = fit_df[,'dir']
  test_dg = test_df[,'drug_cond']
  test_dir = test_df[,'dir']
  data_cols = as.matrix(fit_df[,!(colnames(fit_df)%in%y_colnames)])
  test_cols = as.matrix(test_df[,!(colnames(test_df)%in%y_colnames)])
  model_dg <- cv.glmnet(data_cols,fit_dg)
  model_dir <- cv.glmnet(data_cols,fit_dir)
  pred_dg = predict.cv.glmnet(model_dg,test_cols,s = "lambda.min")
  pred_dir = predict.cv.glmnet(model_dir,test_cols,s = "lambda.min")
  dg_acc= 1-(sum(abs(test_dg-pred_dg))/100)
  dir_acc = 1 - (sum(abs(test_dir-pred_dir))/100)
  running_acc_dg= append(running_acc_dg, dg_acc)
  running_acc_dir= append(running_acc_dir, dir_acc)
  plot(model_dg)
}
return_list = list(mean(running_acc_dg), mean(running_acc_dir))
 return(return_list)
}
### running the function on all our 
all_rats<- list(r1,r2,r3,r4,r5,r6,r7,r8)
j=0
for(this_rat in all_rats){
  j=j+1
  this_rat[is.na(this_rat)] <- 0
  mean_acc<-manual_cv(this_rat)
  print(paste0("mean drug prediction accuracy for rat number ",j," is ",mean_acc[1]))
  print(paste0("mean direction prediction accuracy for rat number ",j," is ",mean_acc[2]))
  
}

## [1] "mean drug prediction accuracy for rat number 1 is 0.979039626237733"
## [1] "mean direction prediction accuracy for rat number 1 is 0.949600724794351"

## [1] "mean drug prediction accuracy for rat number 2 is 0.945991039580656"
## [1] "mean direction prediction accuracy for rat number 2 is 0.948509804974545"

## [1] "mean drug prediction accuracy for rat number 3 is 0.940122702928221"
## [1] "mean direction prediction accuracy for rat number 3 is 0.95149771951286"

## [1] "mean drug prediction accuracy for rat number 4 is 0.948218291046863"
## [1] "mean direction prediction accuracy for rat number 4 is 0.948193156404237"

## [1] "mean drug prediction accuracy for rat number 5 is 0.94067503014628"
## [1] "mean direction prediction accuracy for rat number 5 is 0.946832505854745"

## [1] "mean drug prediction accuracy for rat number 6 is 0.942569572329386"
## [1] "mean direction prediction accuracy for rat number 6 is 0.949935067589102"

## [1] "mean drug prediction accuracy for rat number 7 is 0.939184658223903"
## [1] "mean direction prediction accuracy for rat number 7 is 0.949994824428349"

## [1] "mean drug prediction accuracy for rat number 8 is 0.940983237596321"
## [1] "mean direction prediction accuracy for rat number 8 is 0.951103964746394"

Ignore for now